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I. ANALYTIC SOLUTIONS OP THE HEAT EQUATION SUBJECT TO 
CONVECTIVE AND RADIATIVE BOUNDARY CONDITIONS 

A. INTRODUCTION 

During the 60' s, space technology advanced so much that 
the research of the temperature behaviour of bodies exposed to 
a deep space environment became crucial. In particular, 
transient heat or cooling of solids of different shapes by 
convection and thermal radiation was becoming highly important 
in many engineering applications. An example of these 
applications is the temperature distributions of rocket 
motors. An extensive investigation of the problem has been 
conducted and a lot of literature on the subject was published 
during the 60' s and 70' s. A detailed review of most of these 
papers is not intended here; instead a brief summary of the 
major ones will be given. 

As early as 1962, Fairall, et.al.[6] generated a numerical 
solution for the problem using an explicit finite difference 
scheme; this paper served as pioneer work in the area of the 
research. Later, various finite difference schemes were 
devised to deal with the nonlinear boundary condition. The 
main difficulty in these schemes is the appearance of severe 
oscillations in the determined temperature values for high 
heat flux situations. Von Rosenberg [10] proposed a hybrid of 
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an iterative technique and implicit finite difference schemes 
to deal with the nonlinear boundary condition. On the other 
hand, Crosbie and Viskanta [3,4] transformed the governing 
equations into a nonlinear Volterra integral equation of the 
second kind and applied the method of successive 
approximations to solve the integral equation. Milton and 
Goss [8,9] developed some heuristic stability criteria for 
explicit finite difference schemes with nonlinear boundary 
conditions. It turns out that a very restrictive time step is 
required for numerical stability which may result in requiring 
a prohibitive amount of computer time to calculate the long 
time evolution of the solutions. Williams and Curry [12] 
surveyed several methods for treating the nonlinear boundary 
condition in implicit schemes and compared their accuracy and 
efficiency. 

Nonlinearity is commonplace in natural phenomena. 
Unfortunately, a nonlinear problem often doesn't lend itself 
to a closed form solution. The problem of transient heat- 
conduction in a solid becomes nonlinear when the surface of 
the body is subjected to thermal radiation. When energy 
transfers through the wall of a body, two cases arise: 
convection and thermal radiation. The convective heat 
transfer describes the situation where heat is dissipated 
according to Newton's Law of Cooling, which states that the 
rate at which heat is transferred from the body to a 
surrounding is proportional to the difference in temperature 
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between the body and the environment. The boundary condition 
that describes convection is nonlinear except for the case 
where the heat-transfer coefficient is independent of surface 
temperature, which is technically called forced convection. 
The radiative heat transfer is based on the Stefan-Boltzmann 
Law, which states that the heat flux is proportional to the 
difference between the surface temperature to the fourth 
power and the source temperature. Pure radiation or pure 
convection occur whenever one mode of energy transfer 
predominates over the other. 

It is the purpose of this thesis to consider the one- 
dimensional transient heat conduction problem resulting from 
a combined convective and radiative heat flux with the 
objective of determining the surface temperature fields using 
the numerical methods which are discussed in this study. 
Another purpose of this thesis is to explore the limitations 
and difficulties involved in these schemes. References to the 
work done in similar areas are presented to allow the reader 
further investigation. 

Analytic solutions are derived in one dimension. However, 
the resulting solutions are not in closed form, and thus 
impractical to use. Hence, numerical techniques will be 
studied and employed in the computer in an attempt to find an 
approximate solution. Numerical results, found by 
implementing some of the numerical methods discussed below, 
will be presented and compared. In the conclusion, a 
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numerical scheme is proposed as an alternative to the existing 
methods. It is open to the readers for justification. 

Sections 1(C) and 1(D) describe the derivation of the 
integral representations of the one dimensional transient heat 
conduction problem subjected to a combined convective and 
radiative boundary condition in a rectangular coordinate 
system. Two integral transform methods, namely the Laplace 
transform and the eigenvalue expansion, are presented. 
Observation and comparison are made for the integral equations 
to yield some useful information about the solutions. 

In Chapters II and III, numerical methods for the 
solutions of the nonlinear Volterra integral equations of the 
second kind are described. In particular, the method of 
successive approximations and the Runge-Kutta method are 
outlined in detail. A brief remark is given for their 
advantages and limitations in finding solutions to the 
integral equation. 

Chapter IV describes a numerical method which is directly 
applied to the governing partial differential equation. The 
technique is called the finite difference method. It is 
basically a hybrid of finite difference techniques and an 
iterative scheme proposed. A suggestion is made for the 
improvement of the algorithm. 

In Chapter V, numerical results produced by some of the 
discussed numerical schemes are presented. The implementation 
of various methods gave a practical sense of their advantages 
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and limitations. Graphs and tables are set up in such a way 
that a comparison can be made. 

In the next section, a statement of the problem is given. 
In the statement, the basic assumptions, the governing 
equation and the boundary- initial conditions are included. 

B . STATEMENT OF THE PROBLEM FOR OBTAINING THE SURFACE 
TEMPERATURE 

Considering the one-dimensional, transient, conduction 
heat transfer problem with combined convection and radiation 
at its surface, the following assumptions have been made: 

1. One-dimensional heat transfer to a solid of a finite 
length. 

2. The solid medium is pure, isotropic, homogeneous, and 
opaque to thermal radiation. 

3. All thermodynamic and transport properties are 
independent of temperature. 

4. The solid does not contain any heat sources or sinks. 

5. The fluid is transparent to thermal radiation. 

6. The fluid temperature and the ambient temperature are 
constant. 

The non-dimensional form of the governing partial 
differential equation for the temperature U(x,t) and the 
appropriate initial boundary conditions are 



a 2 u m du 

dx 2 Tt' 



0<X<1, 



t> 0? 



( 1 . 1 ) 
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with initial condition 



U(x,0) = g(x) 



(1.2a) 



and boundary conditions 




(1.2b) 



(ii) dU(l,t) . a3U(lft) » -hU* (1 , t) . 



(1.2c) 



Note: Oi and o 2 can be any real number, except both cannot be 
zero at the same time. a 3 is a non-zero real number, and h is 
a positive real number. 

The next section will deal with solving the partial 
differential equation (1.1) with initial and boundary 
conditions (1.2a-c) by the Laplace transform method and the 
eigenvalue expansion method. As an illustration, two special 
cases with specific values of a x , a 2 , a 3 , and h will be 
considered, and the analytic solutions of these cases at the 
surface will be derived. It will be shown that the surface 
temperature satisfies a singular Volterra integral equation of 
the second kind. At the end of the chapter, we will present 
the solutions and indicate some useful information about the 
integral equations. 
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C. THE LAPLACE TRANSFORM METHOD 



In this section, the Laplace transform of equation (1.1) 
with associated boundary conditions (1.2b,c) is first obtained 
with respect to time. The resulting boundary value problem is 
in terms of the Laplace transform of the required solution. 
Next, the equations are solved for the transformed 
temperature, and the solution of the stated problem can be 
found by taking the inverse Laplace transform of the 
transformed solution. From experience, it can be expected, 
the Laplace inversion is of some difficulty. To simplify the 
situation, specific values of a lf a 2 , a 3 , and h are considered 
so that the inverse process is practical without loss of 
generality. It should be noted that there does exist an 
inverse Laplace transform for other cases of a more general 
nature . 

Now, define the transform of the temperature function, 
U(x,t) , with respect to time as follows 



$£[ U(x, t) ] (s) = f"u(x, t)e' st dt = U(x,s) . (1.3) 

Jo 

After the transformation, the temperature function becomes a 
function not only of x but also of the parameter s. Assuming 
that the derivatives with respect to x pass through the 
transform (differentiation can be accomplished before 
integration) , we have 
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(1.4) 



du(x^t) ] . Q . = (~du(x, t) 0 - s t Hi - = dU(x,s ) 

L 53c J ' ' / o 53c 53c 

9P[ d 2u (*/ fc ) ] ( g) = f"^!iii£^_e' st dt = ^ 2f; ( x > s ) (1.5) 
5x 2 /o 5x 2 5x 2 

The rule for transforming a derivative with respect to time 
can be found using integration by parts. Thus, the Laplace 
transform of the derivatives of U(x,t) with respect to the 
transformed variable t is given by 



S£[ au( *' fc) ] (s) = pdu(x, t) e -st dt = s[/(x, s) - C7(x, 0) . (1.6) 
Ot /o Ot 

Now, applying the Laplace transform to the initial-boundary 
value problem (1.1), (1.2a-c) we remove all time derivatives. 
Holding s fixed, we have the following ordinary differential 
equation in x 



d 2 C 7 j x ; s) - su(x,s) 

dx 2 



-g(x) , 0<x<l 



(1.7) 



with boundary conditions 



a ^diU(0^s)_ - a2 (j( 0 ,s) =0, for x = 0 



(1.8a) 
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(1.8b) 



- U ^ — ~ a 3 U(l,s) = -hS£[U 4 (l,s) ] , for x= 1 

Notice that the initial condition, g(x) , is incorporated in 
the ordinary differential equation. In order to solve (1.7) 
and (1.8a,b), we roust first solve for the general solution of 
the corresponding homogeneous differential equation and a 
particular solution of (1.7) satisfying (1.8a,b). Now, 
consider the general solution of the homogeneous equation for 
(1.7) , 



U hoJ x ' s ) = A e Kx + Be Kx , 



(1.9) 



where 



* 1,2 = 



( 1 . 10 ) 



which are given by the roots of the auxiliary equation 



A 2 - s = 0. 



( 1 . 11 ) 



In the following paragraph we employ the method of variation 
of parameters to solve for a particular solution of (1.7). 
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Let 



U p (x, s) - U J y 1 (x,s) + U 2 x 2 (x,s), (1.12) 

be a particular solution where 1^ ( x , s ) and U 2 (x,s) are any two 
linearly independent solutions of the corresponding 
homogeneous equation. In this case, choose Ui(x,s) = e /sx and 
U 2 (x,s) = e‘ /sx . The object here is to find v 1 (x,s) and v 2 (x,s) 
such that the following equations are satisfied 

e^vi (x, s) + e'^Vj {x, s) = 0, (1.13) 

)[se' /rx x' 1 (x, s) - /se~'' Tx v 2 (x, s) =- g(x) (1.14) 



By Cramer's rule, 



v'i(x, s) 



9_{x) e'^ 
-2fs 



and 



Vz (x, s) 



gr(x) 

2/s 



(1.15) 



(1.16) 



By integrating (1.15) and (1.16), we obtain 



(X/S ) = „ (0 , 

Jo 2 yfs 



(1.17) 
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v 2 U,s) = f’SMe^dz - v 2 (0,s) 

J ° 2y/s 



(1.18) 



Thus, the general solution to (1.7) and (1.8a,b) is 



that is , 



U(x,s) = U h (x,s) + Up (x, s) 



u.iy) 



U(x,s) = Ae v ' 5x + Be' v/5x + (x, s) + e~ y/ * x x 2 (x, s ) , (1.20) 



where A and B are arbitrary constants and ^(XfS), u 2 (x,s) are 
given by (1.17) and (1.18), respectively. To determine A and 
B, boundary conditions (1.8a,b) are used along with the 
following procedure. The derivative of U(x,s) from equation 
(1.20) is found to be 

^ dx ~ ~ A ^ e '^ x ~ Bv/se _v/5x + % /se' /5x v 1 (x, s) + 

e-/ Sx x' 1 (x, s) + e“' /5x V 2 (x,s) - /se~ y/ * x v 2 (x, s) . (1.21) 
Let x = 0. Then (1.15) and (1.16) give 

v£(0,s) + v 2 (0 , s) = 0 . 



(1.18a), (1.20), and (1.21) then imply 
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a x [Ay/s - By/s + v/SVj. (0, S) ~JSV 2 {0,S)] 
- a 2 [A + B +v 1 (0 ,s) + v 2 (0,s)] = 0 



( 1 . 22 ) 



By rearranging the terms, (1.22) becomes 



A(a lV /s - a 2 ) - B(a lV /s + a 2 ) = 

{oCjVs + a 2 )v 2 (0,s) - (a 1 /s - a 2 )v 1 (0,s) . 



(1.23) 



Similarly let x = 1. Then (1.15) and (1.16) give 
e v ' s v£(l,s) + e _v/5 V 2 (1, s) = 0 . 



(1.24) 



Therefore (1.8a), (1.20), and (1.21) then imply 

Ajse'S* ~ B/se~' /s + v'se'^Vi (l, s) - v /se _ ' /s v 2 (1, s) - 
a 3 [Ae v/s + Be _v/5 + e v/5 v 1 (i,s) + e _v/5 v 2 (l, s) ] = 



- hS£[tf 4 (l, t)] . 



(1.25) 



By a similar manipulation of the terms, (1.25) becomes 
Aiy/se'S* - a 3 e v ' 5 ) - B(v'se" v/5 + a 3 e _v/5 ) = 

(v^ie _v/5 + a 3 e _v ' 5 ) v 2 (1, s) - - a 3 e' /5 ) Vj (1, s) - 



h£[U A ( 1, t)] . 



(1.26) 
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Equations (1.23) and (1.26) form a system of two equations in 
the two unknowns A and B. By Cramer's rule, A and B are as 
follows 



_ numl 

den 



(1.27) 



where 



numl - 

((a lV /s-a 2 ) v 1 (0,s) + (a 1 v's+a 2 ) [v 2 (l,s) -v 2 (0, s) ] } {^se'^+a 2 e'^) 
- [ (a lV /s+a 2 )v 1 (l,s) ] (/se' /s -a 3 e v ' 5 ) - h<£. [C7 4 (1, t) ] (a 1 ^s+a 2 ) , 



and 



den = 



(a 2 - a 3 a 1 )/s(e' /5 + e _v/ ®) + (ctjS - a 3 a 2 ) (e^-e"^ 5 ) . 



B = n 'f n2 (1.28) 

den 

where 

num2 = 

{(a lV /s-a 2 ) [v 3 (0, s) -v^l,®)] -(a 1 v^ + a 2 )v 2 (0 / s)}( v ^e' /5 -a 3 e' /5 ) 
- AS£ [U 4 ( 1 , t) ] (a lV ^i-a 2 )+ [ (a 1 v's-a 2 ) v 2 ( 1 , s) ] ( v /se"' / ®+ a 3 e“' /5 ) 

and den is the same. 
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Thus, the general solution of (1.7), (1.8a,b) is given by 
(1.20) where A, B, t!(x,s), and u 2 (x,s) are given by (1.27), 
(1.28), (1.17), and (1.18), respectively. Theoretically, the 
analytic solution of partial differential equation (1.1) with 
initial and boundary conditions (1.2a,b,c) can be obtained by 
taking the Laplace inversion of U(x,s), and thus, the surface 
solution can be found by putting x equal to one in U(x,s) . In 
practice, however, the inverse Laplace transformation process 
is highly unstable in that singularities may exist. Also, the 
transforms are difficult to find. In the next paragraph below 
we consider two special cases where the inversion is feasible. 
In each case, values for the parameters correspond to a 
specific geometrical configuration of a body. 

Case l: a : = 1, a 2 = 0 , a 3 = -1, h = 1 

Thisr set of values corresponds to a "flat plate" with a 
given initial temperature and which is being heated or cooled 
by combined convection and radiation. The term "flat plate" 
is taken here to mean a solid slab of finite thickness which 
is bounded by a pair of vertical lines at ± \ thus having a 
width of 1. Substituting the given values for the parameters 
in (1.20) we obtain the transformed surface temperature 

£7(1, s) = Ae v ' 5 + £e -v ' 5 + (1, s) + e _v ^v 2 (l, s) , (1.29) 
where , 
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[y/gv 1 (0 / s) + v/sv 2 (l,s) -y/sv 2 (0 ( s)] (y/se~^ - e'^) 
v/s(e y/s + e _y/5 ) + sie'/z - e _v/s ) 

+ ~ t/svid/s)] (y/se^ 5 + e^) - A5g [a 4 (1, t) ] Js 
v /s(e' /s + e -y/5 ) + sie^ - e _v/5 ) 



j = [y/iv 1 (0 / s) -y r sv 1 (l / s) - y/gv 2 (0, s) ] (y'se*' 5 + e*' 5 ) 

(a 2 - a 3 a 1 ) v /s(e v/s + e _v/5 ) + (a x s - a 3 a 2 ) (e^ - e" v/s ) 

+ - ASg[C7 4 (l, t)]y/s + [y/sv 2 (l,s)] (y/se - ^ - e -v/5 ) 

(a 2 - a 3 a 1 ) v /s(e v/2 + e -v/5 ) + ( a - a 3 a 2 ) (e^ - e _v/s ) 



v l( l,5) . + V ,(0,s) 

J ° 2 yfS 



and 



v 2 (i,s) - r 1 g(2le^d»'*v,(0,s) 
J 0 2v/s 



(1.32) 



(1.33) 



Now substituting (1.30)-(1.33) into (1.29) and simplifying the 
results gives: 



flg(x') e^’dx' + !\g{x') e'^'dx' 
(e^ + e -v/5 ) + > /s(e v/3 - e _y/s ) 

+ - ASg[tr 4 (l, £)J (e^ t- e-^) 

(e^ + e _v/s ) + v /s(e v/3 - e -y/s ) 
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Suppose the initial temperature is 1, that is, 

g(x) = 1. 



(1.35) 



The boundary conditions associated with the given values of 
the parameters and initial condition (1.35) constitute a 
cooling process. With (1.35), the transformed surface 
temperature becomes 

— (e^-e'^) - h Sl![tf 4 (l, t)] (e^ + e'^) 

17(1, s) = — . (1.36) 

(e v ' 5 + e' v/5 ) + Js[eJs - e" v/5 ) 

If (1.36) is multiplied through by 

1 ^ + e"^ + ys_je^_ - e -v/5 ) j 

Vs (e^ - e-S*) 



and then simplified. 



*7(1, s) 



1 _ &[hU* ll,t) + {e^ + e-S*) 

s y ^( ev/5 ~ e"^) 



(1.37) 



is obtained. Equation (1.37) is ready to be inverted. In 
order to perform the inversion of (1.37), the following two 
Laplace transforms have to be computed 



S2 _1 [ — ] and 

s 



(e^ t e'^) ] 

Vs(e y/s - e" v/5 ) 
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In fact the transforms can be found from any standard Laplace 
transform table. By the convolution theorem, the surface 
temperature in time t is given by 



and by the Poisson summation formula [14], (1.38) can be 

written as 



17(1, t) = 1 - f^ll+2%^ e-* 2n2(t - t) ] [U 4 (1,X) +tf(l,T)] dr . (1.39) 



Hence, the problem of transient cooling of a flat plate by 
combined convection and thermal radiation has been reduced to 
solving a nonlinear Volterra integral equation of the second 
kind . 

Case 2: a x = 0, a 2 = -1, a 3 = 1, h = 1 

This set of values corresponds to the case where a 
spherical body of radius 1 with a given initial temperature 
is being heated or cooled by combined convection and 
radiation. Since the procedures used to solve the problem are 
basically those described in case 1, the mathematical details 
will be omitted and only the main steps will be presented. 
Consider equation (1.20) , the general solution of the boundary 
value problem. The given values for the parameters are first 




[17 4 (1,t) + U{l,x)]ck (1.38) 
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substituted into (1.17), (1.18), (1.27), and (1.28). Then, 

(1.20) is simplified as in the previous case. After a tedious 
calculation, the transformed surface temperature is given by 

f 1 gr(x / ) e'^'dx' - f^g(x') e^'dx' + 

U(l,s) = ^ 12 

(e v/5 _ e ->/8) - y fs(e'/ s + e ~ ' /5 ) 

+ h [ C7 4 ( 1 , t ) ] (e^ 5 + e'^ 5 ) 

(e^ - e + e~v^) 

Suppose the initial temperature is chosen to be 

g(x) = x. 

Boundary conditions associated with the given values of the 
parameters and initial condition (1.41) again constitute a 
cooling process. With (1.41) , the transformed surface 
temperature becomes 

^(i. 5 ) . 1 - Jig[P«(l,t)] (eg - e-^) (1 42) 

s y/s(e + e"*' 5 ) - (e^ - e"' /5 ) 



(1.40) 



(1.41) 



which is now ready to be inverted. In order to perform the 
Laplace inversion of (1.42) the following two inverse Laplace 
transforms need to be computed 



ST 1 !-] and 

s 



rM 

v's(e ' /5 



(e^ - e - *' 5 ) 

+ e -v ^) ~ (e^ 3 



e' v ' 5 ) 
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The first inverse Laplace transform is obvious. However, the 
second one is not so obvious. Details of the derivation of 
the second inverse Laplace transform are given in [1], The 
surface temperature in time t, obtained by inverting (1.42) , 
is 

U{ l,t) =1 - J*[ 3 + 2Y,~ kmi e~ pUt ~ x) ]U 4 (l,T)dx l (1.43) 

where P k is the k— positive root of the transcendental 
equation 

P k = tan p k . (1.44) 

Hence, the problem of transient cooling of a sphere by 
combined convection and thermal radiation has been again 
reduced to solving a nonlinear Volterra integral equation of 
the second kind. As we have mentioned above, one of the 
drawbacks of the Laplace transform method is that there are 
only a few cases in which the transformed solution can be 
practically inverted into the required solution. In the next 
section, the eigenvalue expansion method is introduced as an 
alternative to the above method. One may find the eigenvalue 
method more practical for solving for the analytic solution of 
the heat equation with nonlinear boundary conditions. 

D. THE EIGENVALUE EXPANSION METHOD 

The fundamental idea of the eigenvalue expansion method is 
to transform the given boundary value problem by the 
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eigenfunctions obtained from the associated eigenfunction 
problem. By the completeness theorem (which states that any 
piecewise smooth function can be represented by a generalized 
series of eigenfunctions) we can show that separation of 
variables, i.e., u(x,t) = X(x)T(t), may lead to the solution 
of the problem expressed as an infinite sum of the 
eigenfunctions with appropriate coefficients determined by the 
orthogonality property of eigenfunctions. Applying these 
procedures to the partial differential equation (1.1) and 
initial boundary conditions (1.2a,b,c) yields the following 
main results 

Jl 2 2 fS* }. + p*x{x) - o, 0 < x < 1 (1.45) 
dx 2 



with boundary conditions 

d ^ >) - « 2 *(0) = 0, 

and 



dX( 1) 
dx 



a 3 X(l) = 0 . 



(1.46) 



(1.47) 



Parameters and a 2 can be any real number except they cannot 
be zero at the same time. a 3 is a non-zero real number. 
According to the theory of ordinary differential equations, 
the general solution of (1.45) is 
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X{x) = c 1 cos(Px) + c 2 sin(Px) . 



(1.48) 



Applying boundary conditions (1.46) and (1.47) to equation 
(1.48) gives the following system of equations 



Note that boundary value problem (1.45 - 1.47) is in the 
class of Sturm-Liouville problems for which all eigenvalues 
are real and the eigenfunctions corresponding to different 
eigenvalues are orthogonal. Thus, if the parameters in (1.49) 
and (1.50) are specified, there will exist eigenvalues, P n , 
where n = 1,2,..., and the corresponding eigenfunctions, 

X n (x) , such that the temperature function, U(x,t) , can be 
expanded in a Fourier expansion of the form 



<XiPc 2 = a 2 Ci 

(c 2 P - c 1 a 3 )cosP = (c x P + c 2 a 3 )sinp . 



(1.49) 

(1.50) 






(1.51) 



where the Fourier coefficients, U n (t), are given by 




(1.52) 
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Now, taking the finite Fourier integral transform of the heat 
equation (1.1) with respect to X n (x) gives 



ifJo UlX ' C) X ° lx) 




(1.53) 



Performing integration by parts of the right hand expression 
in equation (1.53) and substituting in (1.52) yields the 
following ordinary differential equation for U n (t) 



With boundary conditions (1.46) and (1.47), the right hand 
side of (1.54) can be simplified. Then, by the integrating 
factor method, the solution of equation (1.54) can be 
obtained. Hence, the resulting integral equation for U(x,t) 
takes the form of (1.51) with U n (t) solved in (1.54). Lastly, 
by putting x - 1, a nonlinear Volterra integral equation of 
the second kind for the surface temperature U ( 1 , t) is 
obtained . 

As in the previous section, the integral equation for the 
surface temperature will be explicitly determined for two 
special cases: the flat plate and the sphere. Details of the 
derivation of the solution will be produced in the case of the 





(1.54) 
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flat plate, but only major results will be given in the case 
of the sphere. 

Case 1: a : = 1, a 2 = 0, a 3 = -1, h = 1 

As mentioned in section 1(C), this set of parameters 
corresponds to the geometrical configuration of a flat plate. 
Substituting the values of a lr a 2 , and a 3 in (1.49) and (1.50), 
c : equals zero, and (1.50) leads to 

cosp a = p fl sinp fl => -jp = tanp,, , (1.55) 

where cos P n * 0. 

So, the family of orthogonal eigenfunctions are 

X a (x) = cos(p B x) , (1.56) 

where n = 1,2,3,..., and {?„)„.!* is the set of distinct 
eigenvalues which are the roots of (1.55) with the property 

O<P1<P2<P3< • • • 



Next, applying the finite Fourier integral transform of the 
heat equation (1.1) yields (1.54) in terms of X n (x). Using 
the boundary conditions 



dUj 0, t) 

dx 



(1.57) 
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- U{ fa — + U(l, t) = - h U*(l, t) , (1.58) 



and the fact that 



*'( 1 ) + *( 1 ) = 0 , 



(1.59) 



*'(0) = 0 , 



(1.60) 



Xnix) 



= ~ PX(*) 



(1.61) 



produces the following ordinary differential equation for 
U n (t) 

d ^ ( . t) + P 2 a U a (t) = - h X n (l) U*(l, t) . (1.62) 



Note that (1.62) is a first order linear ordinary differential 
equation. We find the solution to be 

U a {t) = tf fl (0)e"* 2t - h X a (l) f dx , (1.63) 

J o 



where 



U a (0) = f 1 g(x)X n (x) dx . 
J 0 



(1.64) 
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Thus, with h = 1, the integral equation for U(x,t) takes the 
form 



U(X,C) =^-! { 



[0,(0) «" ,ie - (l) jr„(i) dt)x„(x) 

Jo 



[ 1 X„ ( X ) dX 
J 0 



where U n (0) and X n (x) are defined by (1.64) and (1.56), 
respectively. Lastly, by putting x = 1, the integral solution 
for the surface temperature U(l,t) is determined to be 



where g(x) is the initial condition, and X n (x), and P n are 
defined as above. 

Case 2: a x = 0, a 2 = -1, a 3 = 1, h = l 

In this case, a spherical body is considered. In a 
similar fashion, the family of orthogonal eigenfunctions can 
be found and are given by 




e'^ ot X a [l) f^g {x)X a (x) dx 



- fxtiDe-^uUl,*) dx 

Jo 



”Pl( t*t) rr4 



(1.65) 




X n (x) = sin (P^x) , 



( 1 . 66 ) 
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where n = 1,2,3,..., and P„ is the set of distinct eigenvalues 
that are the roots of 

P fl = tan p a (1.67) 



with the property 

O<PX<P2<P3< 

After applying the finite Fourier integral transform of heat 
equation (1.1) with respect to X n (x) , the following ordinary 
differential equation for U n (t) is obtained 



dU a (t) 

dt 



+ p \u a {t) - - h X n (l)U 4 ( 1, C) . 



( 1 . 68 ) 



Thus, the solution of equation (1.68) is 

U a (t) = - h X a (l) ( 1 e' p2 “ {t ~ x) U 4 (l f x) dr, (1.69) 

Jo 



where 



UJO) = [ 1 g(x)X a (x)dx . 
J 0 



(1.70) 



So therefore, with h = 1, the integral equation for U(x,t) 
takes the form 



- Eli' 



[U„(0)e' p5t - (1) x a (l) dt]x„(x) 

Jo 



Cx 2 n {x) 

Jo 



dx 
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where U n (0) and X n (x) are defined by (1.70) and (1.66), 
respectively. Lastly, by putting x = 1, the integral equation 
for the surface temperature U(l,t) becomes 



where g(x) is the initial condition, and X n (x) and P n are 
defined as above. 

E. REMARKS 

The solution presented above is not complete in the sense 
that the surface temperature is only determined for two cases. 
The solution for other geometrical configurations can be found 
in some of the literature listed in the references, 
specifically 3, 5, 6, and 11. 

The surface temperature solutions which have been derived 
above by both methods fall into the form 

U( 1, t) = 4>(t) - hj*{a + dx , (1.72) 

where F is a nonlinear function of U(l,t) , and c k , b k , a, and 
h are some constants. Equation (1.72) is a nonlinear 




e ( 1 ) J*g (x) X a (x) dx 



f 1 x*(x) dx 
Jo 



- f ^ (l)e' p " <t-t) C7 4 ( 1,t) dx 

Jo 
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Volterra integral equation of the second kind. 4>(t) is a 
function which is usually called the "lag" part of the 
integral equation. The integral in (1.72) is often referred 
to as the "Volterra" part of the integral equation. In 
addition, the piece within the braces of the Volterra part is 
called the "kernel" of the integral equation. As these 
integral equations are being examined, several facts about 
(1.72) are summarized as follows: 

1) . All of these integral equations are singular because as 
x approaches t, the kernel blows up to infinity. 

2) . All of the infinite series satisfy the following 
property : 

If f(t-x) is used to denote an infinite series, 
then lim t _ » f(t) = constant, thus remaining finite. 

3) . The lag part, <j>(t) , and the kernel of the integral 
equations are determined by the geometry of the body 
considered. 

The above "facts" are concluded from the two special cases 
without loss of generality. In each of the next three 
chapters, a different numerical method for solving the problem 
stated in section A will be introduced. Both the method of 
successive approximations and the Runge-Kutta method are 
numerical techniques used to deal with the integral 
representation of the problem, whereas the finite difference 
method is applied directly to the governing equations. 
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II. THE METHOD OF SUCCESSIVE APPROXIMATIONS 



A. INTRODUCTION 

The surface temperature of a body subject to a combined 
convective and radiative boundary condition, as seen in 
Chapter I, is given by the solution of a singular nonlinear 
Volterra integral equation of the second kind. Since the 
integral equation is not in closed form and is nonlinear, 
numerical techniques seem to be the most practical way to 
tackle the problem. Over the past twenty years, a lot of 
research has been done on the numerical solution of an 
integral equation of the form 

*7(l,t) =4>(t) - hf‘{a + dr . (2.1) 

Among the existing numerical methods for solving (2.1), the 
method of successive approximations is the most popular one 
(see [1]). It is based on the idea that the set of successive 
functions defined by 

y ml ( t) =4>(t) - hf t Jc(C-x)F(y n (x)) ch , (2.2) 

where k(t-t) is equal to the term in braces in (2.1), 
converges to a solution of (2.1) in every finite interval of 
time. In the following section, the solution method will be 
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outlined, and at the end of the chapter, general comments will 
be made on the technique. 

B. OUTLINE OF THE METHOD 

Consider the time domain in which integral equation (2.1) 
is to be solved. Suppose the domain is partitioned into N 
intervals. For the first time interval, Oststi , the 
approximate solution of the integral equation can be obtained 
by using the iteration procedure 



until the error between two approximations is less than a 
predefined small number. Next, consider the second time 
interval, t^tstj. In this interval, 



£7 n+1 (l,t) = <| >{t) - Jif t k{t-x) F(U n (l,x) ) dx 

Jo 



(2.3) 



C7 n+1 (l,t) =<J>(t) - hJ‘k(t-x)F(U a (l,x)) dx 



(2.4) 



can be broken into 



U n+ i(l,t) =4 >(t) - hf Cl k(C-x)F(U a (1,t)) dx 

Jo 




(2.5) 



Since U ( 1 , t) is determined for Oststi, the first integral. 



is known (approximately) , and thus the iteration procedure is 
only needed for the second integral. In general, for the i-th 
time interval, the iteration procedure is given by 



U n+ 1 ( 1 , t) =4>(t) - h[ ti ' 1 k(t-x)F(U n (l,x) ) dx 

J 0 

- hf C k( t-x) F(U n (l , x ) ) dx , 



( 2 . 6 ) 



where t^itst*. 

As the procedure continues, the surface temperature is 
found for all desired times. One may notice that as the 
algorithm is carried out, the singularity of the Volterra part 
of the integral equation creates difficulty. Appropriately, 
one has to know the nature of the singularity which the kernel 
possesses'. To illustrate the idea, consider the integral 



This integral is often found in the integral representation of 
the stated problem. Integral (2.7) possesses a singularity 
which can be removed by the use of the transformation 




(2.7) 



Z = a + (b-a) (1-x 2 ) • 



( 2 . 8 ) 
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Then, by using a suitable Gaussian quadrature formula, the 
integral can be evaluated accurately. Normally, one usually 
comes up with an integral with a stronger singularity. 

The iteration procedure outlined above needs a starting 
value. Generally, the algorithm will converge faster to the 
exact solution if the starting value is close to the exact 
solution. Thus, the choice of initial approximation is 
crucial for convergence. Based on the fact that the solution 
of the stated problem is continuous, one can choose the 
temperature at the previous time level as the first 
approximation of the method when a small time step is used. 

The method of successive approximations has been applied 
(in Chapter II) to solve integral equation (2.1). In 
particular, a method used to tackle a simple type of 
singularity, which one may encounter when evaluating the 
Volterra part numerically, has been discussed. Since numerical 
integration is one of the key steps in the method, the choice 
of the numerical integration scheme does affect the overall 
performance of the algorithm. One can improve the accuracy of 
the successive approximations method by appropriately choosing 
a numerical quadrature that can best deal with the singularity 
found in the integral equation. Even though the procedure 
outlined above may seem simple, it has been shown that the 
method is impractical for large times [ 3 ] . 



32 



III. THE RUNGE-KUTTA METHOD 



A. INTRODUCTION 

This chapter considers another way to deal with the 
integral equation for the stated heat conduction problem, 
namely the Runge-Kutta method which was first introduced by 
Crosbie and Viskanta [5] in 1968. The basic idea of the 
method is based on an approximation of the kernel by a 
separable kernel. The integral equation is differentiated 
with respect to time and transformed into a nonlinear 
differential equation. The Runge-Kutta method is a well known 
numerical scheme for solutions of ordinary differential 
equations. In order to employ the method the surface 
temperature at a desired time must be determined. The order 
of approximation of the method is determined by the order of 
the ordinary differential equation. What differentiates this 
method from the other numerical schemes is instead of solving 
an integral equation directly, the Volterra integral equation 
is first reduced to a system of nonlinear ordinary 
differential equations and then solved numerically. The 
method is not exact since the approximation of the kernel is 
not practical if time steps are small. The accuracy of the 
approximation of the kernel increases with time and order. In 
the next section, the method will be outlined in detail as it 
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is applied to the integral equation (1.72). In addition, as 
an example, the formulas for the third and the fifth order 
versions of the method will be presented explicitly. 

B. OUTLINE OP THE METHOD 

Consider the integral equations derived in Chapter I. 
Generally, the integral representation for the dimensionless 
surface temperature, U(l,t), of the body that we have 
considered can be written as 

U( l,t) =<t>(t) - f ft *(t-T)m<l.T)) dx , (3.1) 

Jo 



where 



= Po ♦ EL • 



(3.2) 



The function F(U(l,t)) is the surface heat flux; the a k 's and 
Pfc's are eigenvalues and coefficients, respectively. As shown 
in chapter 1, the infinite series k(t-x) has the property 

lim t _Jc(t) = p 0 . (3.3) 
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This is a necessary condition for an integral equation to 
which the method is applied. Now, assume is a bounded 

differentiable function. The N— -order approximation of k(t-T) 
is given by taking the first N terms of the infinite sum. So, 
(3.2) becomes 

k(C-x) « p 0 ♦ YL • (3 . 4 ) 

Substitute (3.4) in (3.1) and let 

I k (t) = e^J'e^FiUil'X)) dx . (3.5) 



U(l,t) becomes 

*7(1, t) =<|>(t) " P 0 J*F(U(l,x)) dx - PjcJ*(t) • (3.6) 

Differentiating with respect to time, equation (3.6) becomes 

U (1) {l,t) =<J> (1) (t) - p 0 F(I7(l, t) ) 

-EL? k [F(U( 1, t)) - a 2 k I k (t)) . (3.7) 

*7 (2) (l,t) = 4> <2) (t) - p 0 F‘«(^d,t)) 

- JXi P*[F ll) (Z7U, t)) -«|F(tf(l,t)) +a*J*(t)] , (3.7a) 



35 



u {m) { i,t) =4> <jn) (t) - PqF^- 1 ’ (u(i, t) ) 



■ ELi (-D^F^ 1 -^ (U(i, t) ) + (-I) tn al m l k (t)] . (3.8) 



In general, the N— -order approximation of the surface 
temperature U(l,t) is determined by assuming that k(t-T) in 
(3.1) takes the form of (3.4) , in which only the first N terms 
of the infinite sum are considered, so that U(l,t) is 
approximated by (3.6). Then, by performing N+l 

differentiations of (3.6), the resulting system . of 
integrodif ferential equations obtained by substitution of 
1, ..., N+l for m in (3.8) is found to be 

= <J> (1) (t) - P 0 F( C7(l , t) ) 

-ELe k [F(U( 1,0) -alJj^t)] , (3.9) 

U™(l,t) = <|> <2) (C) - PoF* 1 ) (ET(1, t)) 

" T!Li M F<1) (U(l, t) ) - a k F(U(l, t)) + aU k (t)] , (3.10) 
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U m (l,C) = 4> (jN1 ( C) - p 0 F (w_1) (£7(1, t) ) 

- EL mi, t) ) 

* (-l)"afj*(t)] # 



(3.11) 



u {n *d (i, t) =4> ( " +1) (t) - p 0 F ( "> (tfd, t) ) 

- EL P*tEL mi. t) ) 

+ ( -l) t, * 1 a\ ir * 2 I k (C) ] . (3.12) 

To eliminate integrals I k (t) , where k=l,...,N, from (3.12), 
consider the first N derivatives of the surface temperature 
which are given by (3.9)-(3.11) . Rearrangement of the terms 
in (3.9) -(3.11) yields 

£Lcc|p*J k (t) =tf (1) ( l,t) ~ 4> (1) ( t) + P 0 F(tf(l, t)) 

+ EL P Jt F(tr(l, t)) , (3.13) 



"ELi a *M*(t) = tf ( 2 ) (l,t) -<t> ( 2 ) (t) + P 0 F‘«UT(l,t)) 

+ EL p J F<1) {U{1> C) ] - a If(U(1,C))], (3.14) 
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♦ EL P»(EL I'D V 2il F<“-’- J > (£7(1, C) )] , (3.15) 



(-l)"* 1 JLi - V" (1. t) - <I> ,M (C) + PoF'"- 1 ! (£7(1, t) ) 

* EL P*{EL <— (W. t) >] • (3-16) 



In matrix representation we obtain. 



Pi«t 


P 2 «2 






Ji(t) 




B 1 


-P 


””P 2 ®2 


"P*t 




j 2 (t) 




*2 


(-l)"p ia f- 2 


• ••• 
(-l)"P 2 af- 2 ... 


(-D^r 2 




Vi ( 




B W -1 


(-D^Piof 


{-l)" +1 P 2 af ... 


{-D^P^ 




■^w( 







where are defined as 

B 1 = £7 (1) (1, t) + p 0 F(tf(l, t) ) 

♦EL PjtFUTd, t)) , (3.18) 
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b 2 = u< 2) (i, t) - <J> (2) (t) + p 0 f* (1) (u( 1 , t) ) 

+ IXi P*I F<1) (U(1 ' C) ) - a 2 k F(U(l, t) )] , (3.19) 



b n-i = U lN ~ 1] (1/ t) " (t) + Po^' 25 (ff(l, t) ) 

+ ET-i P 4^-o <“ 1 ) (Z7(l, t) )] , (3.20) 

B n =UM (l,t) -<J) (M (t) + PoF^- 15 (17(1, t)) 

+ (-D^F^-^MCTd^))] . (3.21) 

Now, let 



Pl« 2 l 


P2CC2 ••• 


Pw®w 


- Pl«l 


-P 2 «2 


-p«alr 


<-u*Mr 2 


(-l)"p 2 af ’ 2 ... 


(-l)^ 2 " 


(-D w+i p ia f 


(-l)" + 1 p 2 af ... 


(-i)" + 1 p*a; 



By Cramer's rule, I : (t) , . . . , I N (t) can be expressed as a 
quotient of N x n determinants, given by 
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I 2 (t) = 



^N - 1 ( ~ 






B 1 


p2 a 2 




b 2 


-P 2 «2 


-PX 


Bn- i 


• ••• 
(-i)"p 2 af’ 2 ... 


(-D^r 2 


Bn 


(-i)" +1 P 2 af ... 


(-D^P^ 



■* 2 t *2 — 
Pet (A) 



Pi<*i Si P 2 a 2 •• Pw«w 

“Pl<*l -®2 “P3<*3 ••• “Pj/*W 



(-l)"P 1 af"- 2 b n . x (-l)"p 3 af- 2 ... (-l)*P*ajT 2 
(-l)" +1 P ia f b n (-I)^ 1 p 3 af ... (-l)" +1 P*af 





Pet (A) 




Pi«i 


p2^2 ••• -®i 




-Pxtt! 


”P2^2 ••• *®2 


-p*«tf 


(-l)"P ia f- 2 


• ••• • 
(-i)"P 2 af- 2 ... b w _ 3 


(-D^r 2 


(-D^Piaf" 


(-l) w+1 P 2 af ... 


(-l)™P*a 2A 





Pet (A) 




Pl«l 


P 2 0(2 


Pw-i a w-i B 1 


-P^l 


-P 2 «2 •• 


-p*-i«tf-i b 2 


-D^af 2 


(-l)^P 2 af- 2 ... 


(-U^a^ 2 B N . r 


-l)" +1 P ia f 


(-l)" +1 P 2 af ... 


(-1) W+1 P W -1^-1 B w 



Pet (A) 



( 3 . 23 ) 



( 3 . 24 ) 



( 3 . 25 ) 



( 3 . 26 ) 
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where 



Det (A) 



Pl«l 


P 2 «2 




-Pl« 4 l 


- P 2 a 2 


“Pw®^ 


l-DV*” 


(-l)"p 2 af- 2 ... 


(-d^^T 2 


(-i) w+l p ia f 


(-l) Arn p 2 af ... 


(-1) N * 1 P N al t> 



(3.27) 



Next, by using only the fundamental properties of 
determinants, (3.23)-(3.26) can be simplified as a quotient of 
(N-l)x(N-l) determinants (see appendix-A for further detail) . 
Define 



Det {A') = 

-(a\-a\) 

(a 2 -al) 

(-l)"(af 4 -af 4 ) 

( -l ) W+1 (a.™' 2 -a \ N ' 2 ) 

(3.28) 



-(a\-a\) 

(at-ai) 

(-l)"(af- 4 -af- 4 ) 

(-l) N * 1 (al N ' 2 -a 2 1 N ' 2 ) 



-(al-al) 

(al-al) 

(-l) N (al N ' 4 -al N - 4 ) 

( -1 ) " +1 ( al N ' 2 -al N ' 2 ) 
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The formulas are found to be 



Ji(t) = 



{B 2 +CL 2 2 B x ) 


-(al-al) 


... (a£ 0C2) 


(B3-CC2BJ 


ia$-a$) 


( ) 


(B^-t-D^af^Bj 


(-D^(«r-«r) 


C4 N 

1 

r 

rH 

1 


(B^-C-D^af-'BJ 


(-l) w (af 2 -af 2 ) 


... (-I) w+1 (a 7 ' 2 -af' 2 ) 


P x ai Det{A') 



(3.29) 



and for k = 2, N-l, 

(continued on next page) 
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I k (t) 



-1 1 * 



(- 1 ) 



-(a?-aLi) 

(«i-« 4 * 4l ) 

(-1 ) N («? w ’ 4 -aIV) 



-(al-t-aLj) 



( Bj + (X j(. . i Bj ) 
(Bj-aJ.iB,) 



(-l) N *M«r 2 -aM 2 ) 



( -1 ) N ( «*V -«*V ) (B^ - ( -1 ) "aVjB, ) 

( -1) W<1 (al^-al 1 ! i 2 ) (B w - ( -1) ^ai^Bj) 



(a 4 * 2 -a k *i) 



(-l)»(a& 4 -ajft 4 ) 



- (aw -a I*i) 
(a«-a*.i) 



(-D^Mar^-alV) 



P w ajv Bet (/l 7 ) 



(3.30) 
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and 



■Ztf(t) - 




(3.31) 

Thus, I k (t) , where k = 1, N, in equation (3.12) can be 

determined explicitly by formulas (3.29) through (3.31). 
Hence, the integrodif ferential equation (3.12) is reduced to 
a (N+l) — order nonlinear ordinary differential equation 

U (N+1) { l,t) =<t> ( " +1) (t) - P 0 F (W) (£7(1, t) ) 

- EL P*iEL 1) ) 

+ (-l)" +1 ar +2 J*(t)] , (3.32) 

with the initial values 

CT(1, 0) = 4>(0) , (3.33) 

U (1) (1,0) = <J) (1) (0) - P 0 F(Z7(1,0) ) 

- PjtlFUTd^) )] , (3.34) 
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U (2) (1,0) =4> (2 >(0) - PoF^Mird/O) ) 



' EL P*[^ (1) (^(1/0)) - ajF(CT(l,0))] f (3.35) 



U {N) ( 1,0) =d> (M (0) - p o fW-D (^(1,0)) 

" EL ( - 1)i a! i F ( "- 1 - i) ( tr( 1 , 0) ) ) , (3.36) 



which can be obtained by putting t = 0 in (3.9) through 
(3.12), andlj(t), ..., I N (t) are determined by formulas (3.29) 
through (3.31). 

m- 

The Runge-Kutta method is then applied to the nonlinear 
ordinary differential equation (3.32) with initial conditions 
(3.33) through (3.36). Hence, the order approximation of 
the surface temperature will be given by the numerical 
solution of a (N+l)— order nonlinear ordinary differential 
equation. 

As an example, in the next section, formulas for the third 
and the fifth order approximation of the method will be 
presented. Details of the derivation of the equations will 
not be produced, and only major results will be given. 
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C. THE THIRD ORDER APPROXIMATION 



The third order approximation of the surface temperature 
of (3.1) is given by the numerical solution of the following 
fourth order nonlinear ordinary differential equation 

C7 (4) (1, t) = <j) (4) ( t) - (p 0 + Pi + p 2 + P 3 ) F (3) (£7(1, t) ) 

+ (Pi#? + P 2 «2 + P 3 <* 3 ) F {2) (£7(1, t) ) 

- (P ia ! + P 2 a 2 + P 3 a 3 ) F (1) (£7(1, t) ) 

+ (Pi«! + P 2 «! + p 3 af)F(£7(l, t) ) 

- (Pia?! - ! ( t) + P 2 alJ 2 (t) +P 3 a|j 3 (t)) (3.37) 

with initial conditions 

£7(1,0) = <j>(0) , (3.38a) 

C7 (1) (1,0) = <J> (1) (0) - (P 0 + P x + P 2 + P 3 )F(4>(0)) / (3.38b) 

£7 (2) (1,0) = <J) (2) (0) - (p 0 + p 3 + P 2 + P 3 )F (1) (4> (0) ) 

+ (p x ai + P 2 al + P 3 af) F(<J>(0) ) , (3.39) 



and 
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tf (3 >CL,0) = <t> (3) (o) - (p 0 + p, + p 2 + P 3 )f (2) (<t>(o) ) 



ii(t) , 



where 

B. 



+ (P x «? + P 2 <x" + p 3 aa) f (1) (<J> ( 0) ) 

- (Pxai + P 2 a| + P 3 a^)F(<|>(0)) . (3 



I 2 (t), and I 3 (t) in (3.37) are 



I x (t) - - 



(B-j+alBj - (a|-a|) 
(B^atBj (a3-a 2 ) 



p x ai Det(A') 



(3 



J 2 (t) = 



(a? - *!) 

(at -a*) (B 3 -a|s 1 ) 



P 2 a 2 Dec (A') 



(3 



J 3 (t) » - 



(«?-al) (B 2 +alB 1 ) 
(«l-«2) (Bj-atBj 



P 3 a| Pet (A') 



(3 



= tf (1) (1, t) - <J) (1) (t) + (p 0 + P x + p 2 + P 3 )F(P(1, t) ) 



. 40 ) 



• 41) 



.42) 



.43) 
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B 2 = £7 (2) {l,t) - <J> <2) (t) - (P 0 + Pi + p 2 + P 3 ) F (1) (£7(1, t) ) 



+ (Pi«i + P 2 al + P 3 «l)F(ff(l, t) ) , 



B 3 = tf (3) (1, t) - <t> (3) (t) + (p 0 + Pi + P 2 + P 3 ) F (2) (U( 1, t) ) 
- (P x ai + p 2 «i + P 3 a 3 ) F (1) (U( l, t) ) 

+ (PiCC? + p 2 a 2 + p 3 a 3 ) F(J7(l, t) ) , 

and 



Det{A') 



- (a 2 -ai) - (af-ai) 

(a|-al) (a|-at) 



D. THE FIFTH ORDER APPROXIMATION 

The fifth order approximation of the surface temperature 
of (3.1) is given by the numerical solution of the following 
sixth order nonlinear differential equation 
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C7 (6) (l,t) =d > (6) (t) - (P 0 + Pi + p 2 + P 3 + P 4 + P 5 )F< 5 >(£7(1, t)) 

+ (PiCCi + P 2 a 2 + P 3 a 3 + P 4 a 4 + P 5 a|) F (4) {U(l, t) ) 

- (p.a? + P 2 a| + P 3 a 3 + P 4 a$ + p 5 a|) F™ (U(l, t) ) 

+ (Pjaf + p 2 a! + P 3 a 3 + P 4 a 4 + P 5 as) F (2) (U(l, t) ) 

- (Pitt? + P 2 a§ + p 3 a 3 + p 4 a® + P 5 af) F (1) {U{1, t) ) 

+ (PiaJ 0 + p 2 a 2 ° + P 3 a 3 ° + P 4 a 4 ° + p 5 as°) F(U(l, t) ) 

- (p 1 ai 2 J 1 (C) + p 2 a 2 2 J 2 ( t) + P 3 a 3 2 J 3 ( t) 

+ P 4 a$ 2 I 4 (t) +P 5 a| 2 J 5 (t)) , (3.44) 

with initial conditions 

£7(1/0) = <J>(0) / 

£7 (1) (I/O) = 4> (1) (o) - (p 0 + p x + p 2 + p 3 + p 4 + p 5 )F(<M 0 )) , 

ff (a, (i,o) = <J> (2) (o) - (p 0 + p x + p 2 + p 3 + p 4 + P 5 )f (1) (<J>(o) ) 
+ (PiOC 2 + p 2 ai + P 3 a 2 + p 4 a 4 + P 5 a 2 ) F(<|) (0) ) , 
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tf (3) ( 1,0) = 4> (3) (t) - (P 0 + p x + p 2 + p 3 + p 4 + p 5 )F (2) (<j>(o) ) 
+ (Pja 3 + p 2 a? + p 3 a| + P 4 a 4 + P 5 a|) F (1) (<t> (0) ) 

- (P x at + p 2 as + P 3 CC 3 + P 4 a^ + p 5 as) F(<J> (0) ) , 

u {4) (1, 0) = <t>< 4 >(0) - (P 0 + p, + p 2 + p 3 + p 4 + P 5 )f <3) (4>(o) ) 
+ (P x a? + p 2 al + p 3 aa + P 4 a 4 + P 5 al) F <2) (<j> (0) ) 

- (p x at + p 2 d 2 + P 3 a 3 + p 4 at + P 5 a|) F (1) (<}> (0) ) 

+ (Pjaf + p 2 af + P 3 af + p 4 a 4 + P 5 af) F(4> (0) ) , 

£7< 5 ) (1,0) = <j> <5) ( 0 ) - (P 0 + Pi + p 2 + P 3 + P 4 + P S )F <4) (4>(0) ) 

+ (P x ai + p 2 az + P 3 a3 + P 4 a 4 + p 5 al) F (3) (<J> (0) ) 

- (Pi«l + P 2 «2 + P 3 «3 + Mt + P 5 as) F <2) (4>(0) ) 

+ (P 1 a{ + P 2 a| + P 3 a 3 + P 4 a! + p 5 a 6 5 ) F<*> (<J>( 0 ) ) 

- (Pia? + p 2 a| + P 3 a| + P 4 a? + p 5 a|) F(<t> (0) ) . 

lift), I 2 (t), I 3 (t), I<(t), and I 5 (t) in 3.44 are 
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Jl(t) = 



J 2 (t) = 



J 3 (t) 



J 4 (t) 



(B-j+afBJ 


-(af-af) 


- (af-af) 


-(af-af) 


(B 3 -a^B 1 ) 


(af-af) 


(af-af) 


(af-af) 


(B^afBj 


-(af-af) 


- (af-af) 


- (af-af) 


(Bs-alBj 


(«!-«!) 


(af-af) 


(af-af) 




Pjaf Dec (A') 




-(af-af) 


(Bj+afe) 


-(af-af) 


-(af-af) 


(«l - «3) 


(B 3 -a 3 B 1 ) 


(af-af) 


(af-af) 


-(af-af) 


(B 4 +afB 3 ) 


-(af-af) 


- (af-af ) 


(af-af) 


(Bg-afBj) 


(af-af) 


(af-af) 




P 2 af Dec (A') 




-(af-af) 


-(af-af) 


(B 2 +afB 1 ) 


- (af-af) 


(af-af) 


(af-af) 


(B 3 -afB 1 ) 


(af-af) 


-(af-af) 


-(af-af) 


(B 4 +afB 1 ) 


- (af-af) 


(af-af) 


(af-af) 


(B 5 -afB 1 ) 


(af-af) 




P 3 a 3 Dec {A 1 ) 




- (af-af) 


- (af-af) 


-(af-af) 


(B 2 +afB 1 ) 


(af-af) 


(af-af) 


(af-af) 


(B 3 -afB 1 ) 


-(af-af) 


-(af-af) 


-(af-af) 


(B 4 +afB 1 ) 


(«;-«!) 


(af-af) 


(af-af) 


(Bg-afBj 



P 4 af Dec (A') 



( 3 . 45 ) 



( 3 . 46 ) 



( 3 . 47 ) 



( 3 . 48 ) 
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(3.49) 



J 5 (t) 



- (gf-g 4 ) -(gf-gf) -(gf-gf) ( B 2 +alB x ) 
(gf-gf) (gf-gf) (gf-gf) (^-af-Bj 
-(gf-gf) -(gf-gf) -(gf-gf) (S^affi!) 
(gf-gf) (gf-gf) (gf-gf) (-^-gf-Bj 
P 5 a| Det{A') 



where 



Det{A') 



-(gf-gf) 


-(gf-gf) 


-(gf-gf) "(gf-gf) 


(gf-gf) 


(gf-gf) 


(gf-gf) 


(gf-gf) 


-(gf-gf) 


-(gf-gf) 


"(gf-gf) "(gf-gf) 


(gf-gf) 


(gf-gf) 


(gf-gf) 


(gf-gf) 



B 1 = V (1) (1, t) - <J) (1) ( t) + (p 0 + Pi + P 2 + P 3 + P 4 + V 5 )F(U(1, t) ) 

B 2 = U< 2 1(1, t) - <|> (2) (t) + (p 0 + p, + P 2 + p 3 + p 4 + P 5 ) F (1) (£7(1, t) 

- (P x gi + P 2 *2 + P 3 a 3 + P 4 af + p 5 gi) F(U{ 1 , t) ) , 



b 3 = £7< 3 > (1, t) - <i>< 3) (t) + (p 0 + p, + p 2 + p 3 + P 4 + P 5 ) F (2) (17(1, t) 
- (Pi«f + P 2 «l + P 3 «! + P 4 ®4 + Psgf)^ 1 ’ (£7(1, t)) 

+ (P x gf + P 2 g 2 + p 3 g 3 + P 4 g 4 + p 5 g|) F(U(1, t ) ) , 
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